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Abstract 

As a first result, an explicit analytic expression for the density of the sum of two T 
random variables is derived. Then a computationally efficient algorithm to numerically 
calculate the linear combination of an even number of Xr random variables is developed. 
An explicit expression for the error bound is obtained. As an application, the exponential 
!— 1 dependence of the eigenvalue fluctuation probability of a random matrix model for Ad 

supergravity with N scalar fields is found to be of the asymptotic form exp(— 0.357V). 

> 

1. Introduction and Review 

1.1. The Sum of Chi Square Random Variables 

od 

Many fields of research profit from progress in random matrix theory. Especially relevant in 
string theory, random matrix theory has proven an invaluable tool for making statistical 
claims about the existence of stable vacua. We will revisit one of these applications in 
# ^ section 5. 

^ Linear combinations of x 2 random variables occur frequently in the study of random matrix 

models, such as the spectral analysis of the Wishart ensemble. For a full understanding of 
a random variable, its distribution is of course of utmost importance. The great number of 
related publications proves the ongoing interest in this topic. 
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1.2. The Existing Literature 

As already mentioned, there exist many approaches to calculating or approximating the 
distribution of a linear combination of x 2 random variables. For our purposes, one exclusion 
criterion is that the method has to be quantitatively accurate, which means that we need 
an explicit expression for the error. This excludes fits to other distributions, as for examples 
the three moment fit suggested by [Solomon and Stephens, 1977] or the more recent mixture 
approximation by [Bruce G. Lindsay and Basak, 2000]. Let us discuss the techniques which 
pass this requirement. 

One brute-force method is of course to draw a couple of million random variables, multiply 
them with the weights and calculate the distribution empirically. This works well when 
one is interested in the quantiles that lie around the bulk of the mass. In our example, 
however, the very low quantiles are of interest, and we simply do not have enough statistics 
to calculate those. Note that importance sampling can — unfortunately — not be used, since 
we do not know the complete distribution. 

Another often-used ansatz is to expand the moment generating function in a series 
representation and inverting term by term. One such example (for gamma random variables) 
is given in [Moschopoulos, 1985]. Unfortunately, the prefactors in the expansion cannot 
be calculated in polynomial time. Another example (for chi squares) is [Moschopoulos, 
1984] where, again, the prefactors are hard to compute. A more recent example is [Castaho- 
Martinez and Lopez-Blazquez, 2005], which uses Laguerre polynomials. The same problem 
persists here; in addition, the error bound given suggests that obtaining an accurate result 
will require a very large number of terms. 

As a third group, there are techniques that use numerical integration to obtain the density 
function, see for example [Fleiss, 1971] and [Akkouchi, 2005]. Unfortunately, in both papers, 
the number of dimensions over which we have to integrate for n terms is n — 1. While this 
might be feasible for small n, it is not for our application — just sampling two points in 
any direction grows as 2 n , and even using sophisticated Monte Carlo or Gibbs sampling 
algorithms does not help when n is of order 1000. Another method, which has only one 
integral to calculate, [Davies, 1980], has the problem of a highly oscillatory integrand that 
has to be integrated over R, which seems impossible to do to the required accuracy. 

To summarize, none of the above-mentioned techniques is feasible in our regime of 
interest — either due to computational restrictions or because the answer is not precise 
enough. In this paper, I develop a method which is both fast, arbitrarily accurate and 
computationally efficient, especially in the low quantile regime that we are interested in. 

2. Analytic Results 

2.1. The Linear Combination of Two Gamma and x 2 Random Variables 

Definition 1. Let X ~ T(a,/3) be a gamma- distributed random variable with shape pa- 
rameter a £ R>o and rate parameter (3 £ R>o- Then its density function fx is given 
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by 

fx(x)=e(x)J^-x a - 1 e-^ . (1) 
r(a) 

Note the following 

Remark 1. // X ~ r(a, /?), then cX ~ c/3) for c £ R >0 . 

This is why /3 is also called inverse shape parameter. This immediately leads me to the 
first 

Theorem 1. Let Z = X + Y where X ~ r(ai,/?i) and Y ~ r(a 2 ,/3 2 ) ore independent 
Gamma distributions. Then the density function fy is given by 

fz(z) = e(z) ^ P 2 , z^ + ^~ 1 e~^ l ¥ l {a l -a 1 + a 2 ; (ft - fa)z) , (2) 
T(ai + a 2 ) 

where iF\ is a confluent hypergeometric function (Kummers function of the first kind). 
Proof. We have to calculate the convolution fz = fx* fy over R. 

r(aijr(a2j Jr 

= /^rg 2 e -fc*g( z ) r _ x )« a -ie(A-ft)*ds , 

r(ai)r(a 2 ) Jo 

which can be further simplified by substituting x = zt dx = zdt, so 
(fx*fy)(z) 



r(a 1 )r(a 2 ) 



2 a ~foz 



r 1 
Jo 



r(ai)r(a 2 ) r(ai + a 2 J 

= 9(z) Pl P t 2 ^ +a "- 1 e- fe2 iF 1 ( ai ; ax + q 2 ; (ft - ft)z) , 
r(ai+a 2 ) 

where we have used the integral representation of jFi in the second to last line. □ 

Observe that equation 2 is indeed symmetric under exchange of the subscripts i o 2 , as 
expected. To the best of the author's knowledge, this result is not available in the literature. 

Definition 2. Let X ~ \\ ^ e a chi-square random variable with k degrees of freedom. Then 
the density function fx is given by 

k ^ x 

f x (x) = 9(x fl 6 * . (3) 

2*r(f) 
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Remark 2. If X ~ x£, then aX ~ T (J, /or a G R >0 . 

As a second important result, which will find its application later, we conclude the 
following 

Corollary 1. Let X,Y ~ xl ^ wo chi-square random variables with k degrees of freedom. 
Let Z := aX + bY , a,b £ R>o- Then the density function fz is given by 

i_ k -p ( 1 _ k\ 

fz{z) = 6{z) TAaW [-Sab J 6 V§ l^T^ ' ^ 

where \ v is the vth order modified Bessel function of the first kind. 

Proof. This follows from theorem 1, remark 2 and from the identity (known as Rummer's 
second transform) 



lFl (a; 2a; x) ^ e§ (f ) Q + a) I a _ § (|) . □ 

Equation 4 struck my interest for the following reason: for small z, given a and 6 are 
of the same order, the prefactor in the modified Bessel function is small, so a promising 
ansatz for an approximation should be a Taylor- Maclaurin series for \ v . The exponential 
then restricts most of the mass to lower values of z, which, in addition, suppresses the error 
introduced from the expansion. Furthermore, the leading nontrivial expansion term of Io(-z) 
is of order 0(z 2 ) and for general \ v only every second order monomial is present. All this 
will be quantified in the next chapters. 

At first, though, we need some more results that will come into play later. 

2.2. Convolution Algebra 

Definition 3. Let * : C(R) x C(R) ->• C(R) be defined as 



if * g){z) = 9(x)9(z - x)f(x)g(z - x)dx = f(x)g(z - x)dx 
Jr Jo 

for f,geC(R). 

This map is well-defined, as follows from [Stein and Weiss, 1971], Th. 1.3, and could 
well be extended to spaces such as L 1 (1R), but we do not need more here. Important is the 
following 

Theorem 2. Let P N := {e biX x ni : hi G R, rij £ M ,nj < N} where JVel. Let further 
f, g G Pat: f(x) = e ax x n , g{x) = e bx x m and define 



(/ * g)(x) ■-- 



: a = b 

(/ * g){x) :a^b 



Then ((Pn)^,*) forms an algebra over R ViV G M, where * extends canonically to the 
linear hull. 
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Proof. The case a = b is trivial. Analogously to the proof of theorem 1 , we have 



(/ * g)(x) = 9{x)e bx T{ ™. + ^ )r( ^ ± l) ^+" +1 1 F 1 (m + 1; n + m + 2; (a - 6)x) . 
1 (n + m + 2) 

Using that, for integers r,s£ IN>o, r < s, we have (see [Wolfram Web Resources]) 

[to k ^~^ to kl(2-s) k ) ' 

we note that the lowest occurring power of z is of order 0(z 1 ~ s ) or, including the previous 
result, a constant. The highest occuring power is of order 0{z~ r ) resp. 0(z n ) or 0(z r ~ s ) 
resp. 0{z m ). □ 

Note that the use of * over * is just a technicality. We will see that we can use * as long 
as our exponentials are distinct. 

Theorem 2 is an important result, as it claims that convolving two terms does not result 
in higher order monomials in our expression. But we can make even stronger claims. 

Remark 3. Let f\, .../„ £ ~Pn> gi, ■■■,g m E Pftf, all fi,gj having pairwise distinct exponen- 
tials, and ai, a n , b%, b m £ R. Then 31 < (n + m)N, h±, ...,hi £ Pn pairwise distinct 
and ci, q E R: 

i m \ I 

K i=l j=l J i=l 

Proof. Here, the same notation as in the proof for theorem 2 is used. Counting the number 
of independent exponentials l e , we obtain l e = n + m. The maximum number of monomials 
multiplied with an exponential is then max{s — r — 1, r — 1} = max{n, m} < N. □ 

This by itself can still grow very quickly for more than two convolutions, but the next 
result puts a stricter bound on the number of terms if we convolve multiple times. 

Theorem 3. Let fi, f q £ (Pjv)r so that fi has ni distinct exponentials, and a±, a q £ R. 
Then 31 < N J2l=i n «> h\, ...,hi £Pn pairwise independent and ci, q £ R: 

l 

(ai/i * ... * a q f q )(x) = ^2ahi(x) . 

i=l 

Proof. Let us consider the first convolution. We first note that the number of distinct 
exponentials l e in remark 3 is l e = n\ + 712, even if we replace the monomials with arbitrary 
polynomials. The rest then follows from induction and the fact that ((Pjv)ir > *) — being an 
algebra — is closed. □ 

With these basic facts, let us now discuss the actual algorithm to calculate the distribution 
of a linear combination of chi square random variables. 
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3. The Density Function of a Linear Combination of x 2 Random 
Variables 

3.1. Density Function 

I will first claim the main result and discuss the details in due course. 

As stated in the first chapter, the problem is to calculate the density function of Z = 
Y2=l a iXi, where Xi ~ Xr are hd random variables for some r > and G H>o- 

Let T n : C n — > R[ra] denote the Taylor-Maclaurin map which assigns to a function its 
Taylor polynomial to order n around 0, and R n : C n — > C n its associated remainder map. 

Theorem 4. Let r G 2Mo + 1, Z = J27=i i a i^i + ^jXi), Xi,Yj ~ Xr iid random variables, 
ai,bi G R>o for some n > and a\ < b\ < ... < a n < b n . Let mi, ...,m n G Mo and 
rh := maxjmi, ...,m n }. Then the density function fz is given by 

f z (x) = Ce(x)(T(x) + R(x)) , 

where ^ =1 fi ■= fi*-*fn, 

t W .(*.-^(t w , h )(^-))w 



and 



r 



r / "» 



bi - Qi \ 2 2 



and 

< R{x) < R(x) - T(X) 

*«- 5 *-(T„I | ., + R < )(^.))(x)-T(,) 

i=l 

Finally, \/x < y G R>o, Ri is bounded by 



Ri(x) := R mi Ir_i f ^ — z^aA 
2 2 V 4ai6i / 



~~ (m» + l)! V 4aj6j / §-§ \ 4aA 

Proof. The main statements all follow from corollary 1 and Taylor's theorem, where 
R_i G C°°(R) Vr G 2Mo + 1. The estimate for the remainder is a simple uniform estimate, 

2 2 

using the fact that I^ k \x) is strictly monotonically increasing \/x G R>o,^ > and 
k G Mq. □ 
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It might seem like a restriction that r has to be odd. Note though that the case r even is 
simple, since — as already mentioned — we would only have to convolve terms from P^r (see 
definition 2). Also observe that it is not of great significance that this expansion holds only 
for an even number of terms, as convolving numerically once is easily accomplished. 

In the light of theorem 3 and if r G 21No + 1, the claim is that / is very efficient to 
calculate to high precision, especially if the argument is not too big. To quantify this a bit 
more, we have to relate the magnitude of the error to the prefactors ctj, hi and the abscissa 
x. 

3.2. Error Estimate and Complexity Class 

Let us consider the general case. We are interested in all x G [0, x max ) =: Q, for some 
s max > 0. Let /, g G L 1 ^) n C°°(Q) positive and T n g f g. Then 



\(f*(g-T n g))(x)\= / f(x-y)(g-T n g)(y)dy 

Jn 

< I \\f\\%{9 -T n g)(y)dy = x m&x \\f\\n f (g-T n g)(y)dy . 
Jn Jn 

This of course also holds when we have an additional weight function r\ in the product, 
which we will denote by *„. Assume now we want to have the error ratio below some fixed 
-Rmax > 0. Then one can estimate 

(/ * v (g-T n g))(x) x msK \\rif\\^ S a ri(y){g-T n g)(y)dy 



(/ *r? g)(x) x max \\r]f\\^ f n r](y)g(y)dy 

_ 1 / n ^)(T rag )(y)dy _ 1 |MTng)||? 



In^{y)9{y)^y llwlli 



n 



or 



||»7(T n 5)||?>(l-i2max)|rf • (5) 



For the numerical algorithm, this is basically what one uses to determine how many Taylor 
terms have to be kept to satisfy the inequality. This integration can be performed very fast 
numerically, since the integrand is sufficiently regular and the integration region compact. 

Since this result does not depend on / — which is one of the reasons why the error is 
exaggerated gravely, more of this later — we can inductively determine how many terms we 
have to keep: in the most naive way possible, for m functions, we remain below the error 
bound -R m ax if each convolution separately satisfies the error bound ^rff • 

Just to get the order of the complexity class right, let us consider our special case 

nix) = e~ x and q(x) = It i ( t~x\ for a, b G R>n 

' w yy ' 2-2 \ Aab ) 



Then 

\n . w^\\n\\rji 

loo 



IWT„ 5 )||?> NlfHT^I 
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using Holder's inequality, and applying the same to the right hand side of equation 5, we 
obtain the relation 

\\T n g\\^>(l-R max )\\g\\^ . 
Since \ v is monotoneously increasing, this is equivalent to 

(T„#)(x max ) > (1 - Rma,*)g(x max ) 

or 

(^max) 5; -Rmax<?(-£max) 

Let now y := %max- Using the uniform estimate from theorem 4, 

(n + ly. 2~2 2 2 

For r = 1, e.g. if all our random variables are Xi-distributed, we know that Ig (a;) < 
Io(s) Vx > 0. We need 



, _ e" r log ( — ] < -- <= n > rW" 1 f — 9 

n! \ n I n \x z 



where, in the second step, the inequality n\ > na is used. For large n, it could be replaced 
by n! > n rn for some v — v n G j"^, 1^ . W is the Lambert W function. With this result, we 
finally require 

< ^Rmax =: e~ rmax or roughly n G 0(r max yf ' 



(n + 1)! 



max 7 



Note that these final estimates tremendously overestimate the error and should not be used 
to actually determine the number of terms that have to be kept. Equation 5 gives a far 
more reliable result. 

Combining this with our results from examining the convolution algebra, theorem 3, we 
get the following 

Theorem 5. Let the setup be as in theorem 4 ind r = 1. Then the computational 
complexity to calculate a value x of fz with relative error smaller than it! max £ (0, 1) is 
0(n 2 y 2 (logn-log(i? max ))) ; where y := maxj =1 ,.,r.{^}. 

Proof. Consider the ith convolution. A straightforward algorithm has to convolve Ni ■ i 
terms with N2 • 1 terms. Using the notation from above, an upper bound surely is iVi = N2 G 
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O ^ log( ^™j* ) ^) — 0((log(n)— log(i? max ))y), thus — in a straightforward implementation — 

we will have to convolve the terms pairwise, i.e. i ((log(ra) — log(-R max ))y 2 ) 2 times. Summing 
all steps up, one concludes that the number of terms that have to be convolved are 

/n-l 

O (J2 i ((log(n) - log(i? max ))y 2 



2 



vi=l 



= of (log(n) - log( J R max ))y 4 • ^n(n - 

ylog(/? max ) y 

This is the second main result. In words this means that the algorithm grows at most 
polynomial in the number of terms in our sum, polynomial in y as given above and 
logarithmically in the precision we want to obtain. 

This of course neglects implementation details — I assume for example that the integration 
estimate equation 5 is constant for all terms. The convolution itself can be done by 
multiplying matrices and then sorting vectors, which is more likely to be in O(nlogn), but 
this would only introduce a minor change in our end result. 

Moreover, one finds that for our application (see section 5), the overall growth is rather 
linear in n, so it remains an open question to obtain a better estimate on the complexity of 
the algorithm and to generalize the result to arbitrary r. 



3.3. A Better Error Estimate 

Instead of discussing implementation details, let us look at the error bound. Note first the 
following 

Remark 4. Let £1 C R be a closed interval, f G C°°(£l) : f^ k \x) > Vi £ SI and 
m > n G Mo arbitrary. Then 

(R n /)(x) < ((T m -T„)/)(x) + M m Vrreft , 

where (R m /)(x) < M m Vx G 0. 

Proof. This is a direct consequence of the Lagrange remainder formula and probably 
obvious. □ 

The application for this remark is a bit more subtle. Imagine we want to calculate the 
distribution for n terms. By using the direct Lagrange estimate as stated in theorem 4, 
we will overestimate the error by quite a bit. If this effect is significant, it could lead us 
to expansions to much higher order for both the estimate and the remainder (note that 
the polynomial order of the remainder is always one more than the order of the estimate 
itself). However, given our new error bound, we save ourselves half of the work, namely the 
expansion of the estimate polynomial. I have found this improvement to be very effective. 
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0.5 1 1.5 2 2.5 3 3.5 4 0.5 1 1.5 2 2.5 3 3.5 4 

(c) Absolute errors, few expansion terms (d) Absolute errors, many expansion terms 

Figure 1: The distribution for 20 x 2 random variables. Figure (a) shows the density 
compared to an empirical histogram distribution. Figure (c) shows the absolute 
predicted error vs. the true error, compared to the histogram distribution, (b) 
and (d) is are the same, with more expansion terms. Note that the uncertainty in 
the histogram distribution is — due to computational restrictions — huge. 

4. Numerical Analysis 

Let us now turn to a numerical analysis of my proposed technique. I implemented it using 
Mathematica. As already mentioned, equation 5 is used to decide how many terms to 
keep. The actual convolution in theorem 4 can be implemented very efficiently as a simple 
algebraic method, since we never leave our * algebra. 

4.1. Comparisons with other Methods 

One way of checking our results quantitatively is comparing them to Monte Carlo simulations, 
namely building a histogram distributions from a table of drawn values. This will only work 
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(a) Laguerre expansion 



(b) Numerical integration 



Figure 2: The cumulative distribution for 6 x 2 random variables. Figure (a) shows an 
expansion in Laguerre polynomials as proposed by [Castaho-Martinez and Lopez- 
Blazquez, 2005], with 1, 6, 11 (blue) and 16 (orange) expansion terms. Figure (b) 
shows a numerical integration proposed by [Akkouchi, 2005]. The black line is my 
method, which is exact at the shown resolution. 

in the regime with a lot of mass in the distributions, the tails will not receive enough data 
points to be quantitatively correct. 

In figure 1, you can see the density for n = 20, where the weights were taken from 
equation 7. Depending on the threshold that is taken for the error bound, equation 5, the 
expansion comes arbitrarily close to the empirical distribution. Note that the error bound 
shrinks accordingly and — as expected — we only have an error in the positive y direction. 
Also observe that our error bound still vastly overestimates the true error. 

In figure 2, my method is compared to two more recently published techniques, where 
the six weights are again taken from equation 7. The series expansion (a) can be made 
quite accurate (although the error bound — not shown — is very bad and suggests we expand 
to much higher order), at the cost of exponentially growing computation time. For more 
than 10 terms, it seems impossible to expand to high enough order in any sensible time. 

On the other hand, the numerical integration — performed with Mathematical adaptive 
quasi Monte Carlo method with about 10 5 integrand evaluations — suffers from severe 
convergence problems. 

Just as a comparison: The Laguerre expansion takes about 10 seconds for 16 terms, the 
numerical integration about 250 seconds (for all the values shown — for a single value it is 
about 10 seconds), whereas my method needs about 0.2 seconds. While not representative, 
this should give a rough idea of what to expect; they were all run on one i7-2670QM core. 
While it is true that implementing the Laguerre expansion or the numerical integration in C 
could be orders of magnitudes faster, mind that our method is executed by Mathematica, 
too. 
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5. Counting String Vacua 



5.1. Background and Motivation 

String theory admits some 10 500 vacuum solutions — for an interesting digression, see 
[Susskind, 2003]. They emerge from the parameter space describing the internal compact 
manifold, which string theory requires in addition to our 3+1 dimensional spacetime. This 
parameter space is called the moduli space of super symmetric vacua. At low energies, these 
moduli appear as massles scalar fields with their own equations of motion and potential. 

The question of how many of those configurations correspond to physically stable solu- 
tions is an active field of research. One vacuum configuration — a deSitter spacetime — is 
particularly physically motivated, since it agrees with the recent discovery of an accelerated 
expanding universe. 

In the study of J\f = 1 supergravities, it was shown in [Marsh et al., 2012] that the 
fraction of metastable deSitter vacua shrinks exponentially with the number of scalar fields 
N S> 1. By constructing a random matrix model for the Hessian matrix H of the scalar 
potential and calculating its eigenvalue spectrum, local minima can be found. 

The suggested functional dependence for N fields for the probability of the smallest 
eigenvalue fluctuating to a positive value — which renders H positive definite — was 

P N cx exp (-cN p ) where c,t>0 . (6) 

The Hessian matrix can be block-diagonalized to H = diag(m_,m + ), where m± denotes 
a mass matrix. An important result is the explicit expression for the smallest eigenvalue of 
the mass matrix, 

m 2 = F 2 T± - F 2 S . 



7± and S are given by 
2 



T ± = --u 2 + Kt l K lle -K im ± 



3 



N ITT |2 
Ullb'\ 

b' =2 A b> 



where 



s/3 

to = —\W\ and F oc 



F 1 1 VN 

Here, K denotes the Kahler potential and W the superpotential, which are taken to be 
random functions. U and F are derivatives of W. As shown, the statistical properties of 
T± can be neglected for large and replaced by its expectation value, so T- ~ 1. The 
last term thus gives the leading contribution to the fluctuation probability, i.e. P(<S < 1). 
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0.5 1 1.5 2 2.5 3 3.5 4 10 20 30 40 

(a) (b) 



Figure 3: (a) The Marcenko-Pastur law. (b) Weights for N = 24 and ./V = 42, as determined 
by equation 8. 

The Uny are N — 1 independent and identically normally distributed random variables 
~ A/"(0, 1), while the denominators A 2 / are determined by the Marcenko-Pastur law 

M X ) = 27: XNa^ X{4Na2 ~ X) Where a= VN ' (7) 
It can be justified to take the weights A^ 2 such that 

XI = (XI) : ff® / MP (A)dA = b - . (8) 

This corresponds to numerically inverting the function /mp on the range (0,4), which is 
well defined (unfortunately, this cannot be done analytically), /mp and the weights for 
specific N can be seen in figure 3. 

In order to calculate the fluctuation probability P(<S < 1), we need to know the distribution 
of S, which is a linear combination of iid Xi random variables (chi-square with one degree 
of freedom). To my knowledge, a closed, simple expression for this is not known. In [Marsh 
et al., 2012], a three moment fit ([Solomon and Stephens, 1977]) is used to approximate the 
distribution with a single Xp variable — the quality of which we will discuss later. By then 
fitting the iV-dependence to equation 6, the parameters obtained were 

p = 1.28 ±0.10 and c = 0.29 ± 0.06 . 
In contrast, by simulating the full mass matrix numerically, the values 

p 0.24 and c w 23 
were computed. The difference is significant. 
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Figure 4: Figure (a) shows the N dependence of Pn- Figure (b) shows that the result is 
coherent with the simulation of the mass matrix, as done by [Marsh et al., 2012]. 
Figure (c) shows the computation time needed to calculate the values and (d) the 
relative errors, as determined by the algorithm. 
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5.2. Simulation 



Let us return to the issue of finding the N dependence of the fluctuation probability of the 
mass matrix in our Ad random supergravity. We want to calculate 

P N :=-p(F 2 J2^^<lj ^.V^aiXiKlJ =:F(Y N <1) , (9) 

where X x ~ x 2 (l) and a, := (JV(Agr))- 1 . 

Using the method discussed in section 2, we calculate the pdf fy N for N from 2 to 2250 
and integrate it over (0, 1), which gives -P/v- All the while the simulation keeps track of the 
relative error and ensures that it stays below 0.05. The results can be seen in figure 4. 

The first 400 values were calculated on a single i7-2670QM, while the rest were computed 
on a small cluster with slower individual CPUs (i7-860) and varying load balance, which is 
where the big amount of noise for the computation time and the relative errors originates. 

The simulation agrees very well with the numerical simulations of the mass matrix, as 
can be seen in the low TV" regime figure 4(b). For larger N, figure 4(a), the overall shape 
then follows more and more that of a simple exponential, i.e. tx e~ cN for some c. 

5.3. Data Analysis 

Indeed, by fitting the suggested functional dependence e~ cNP , the residuals show a strong 
characteristic shape, see figure 5. While the overall shape matches pretty well, another 
function, 

e log(a+iV)- C -dAr = (a + iV)e - c - ( iJV ^ (1Q) 

matches the shape a lot better, as can be seen by looking at the corresponding residuals. 
While figure 5(d) still shows a correlation for the residuals, they are now of the same order 
as the proposed error, figure 4(d). This means that this model is correct, within the given 
precision. 

The fit parameters and standard deviations for both models can be found in table 1. 
This is the major result of my application. In brief, it means that the probability of 
finding a metastable deSitter vacuum follows the asymptotic form 




6. Conclusions and Future Work 

The first main result is the analytic expression for the density of a sum of to T random 
variables (theorem 1). With this expression as a building block, I have derived an efficient 
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(c) Residuals for (a) (d) Residuals for (b) 

Figure 5: Figure (a) and (b) show plots of fits (orange) to our data (blue) for both proposed 
models. Figure (c) and (d) show the corresponding residuals. 



(a) Parameters for —cN p 


(b) Parameters for log(a + N) — c — dN 


Estimate Standard Error 


Estimate 


Standard Error 


c 0.31395 0.00011 


a 1.88 


0.06 


p 1.01171 0.00005 


c 0.75652 


0.00011 




d 0.34657360 


4 • 10~ 8 



Table 1: Fit parameters for both proposed models. 
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algorithm (theorem 2) for the calculation of a linear combination of an even number of Xr 
random variables, which, in the case r = 1, is proven to be computationally efficient. 

It remains an open question to refine and extend theorem 3 to arbitrary degrees of 
freedom. Nonetheless, I am confident that the algorithm will stay polynomial in n. 

With these methods, a reliable estimate for the probability of a metastable deSitter vaccum 
in a M = 1 Ad supergravity has been shown to be of the asymptotic form Ne~°' 35N , where 
N denotes the number of scalar fields in the theory. As discussed, this last result gives a 
much weaker bound than the one suggested by [Marsh et al., 2012] or [Chen et al, 2012]. 

This result is significant. One important implication for the existence of flux vacua is the 
following. If we assume that there are M £ O(10 500 ) local extrema, of which M + <C M are 
metastable deSitter vacua, we can conclude that 



log(a + N) - c - dN ~ log N - dN ~ log 



M + 
~M 



Consequently, this description of a 4d random supergravity is plausible — i.e. M + > 1, 
which means that there is at least one metastable vacuum — if —500 log 10 < log N — dN or 
N < 3300. 
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A. More Exact Formulae for x 2 Random Variables 

There are a few more analytic expressions for the pdf / of a sum of x 2 random variables. 
Remark 5. Let X,Y, Z ~ Xi random variables and a,b € R>o- Then 



faX+aY+bz{x) = 0(x)\ — e 2f > F 





jb-a 


na 2 b\ 


ab 




k-a 


4a 2 b\l 


ab 




d ( x )\ tz^i\ H^r e 2aerfi 



where F is the Dawson integral and erfi the imaginary error function. 
Remark 6. Let X, Y, Z,W ~ Xi random variables and a, b S R>o. Then 

faX + aV +a Z +b w{x) = ^)\^-^ X (*> (^*) + h 
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